Current density imaging on the epicardial surface of the heart 
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1 Introduction 

Magnetocardiographic (MCG) mapping has provided 
promising results in locating various intracardiac 
sources, such as the sites of origin of life-threatening 
arrhythmias (e.g., [1]). Magnetocardiographic inverse 
studies with point-like source models, however, are 
generally applicable only for sources confined in a 
small volume of tissue, e.g., at the onset of the ven¬ 
tricular depolarization. After the onset, the electrical 
activation spreads in the myocardium as a wavefront, 
and the use of point-like source models can no longer 
be physiologically justified. Therefore, distributed 
source models provide a more general description of 
the current sources in the heart. In this study, the ap¬ 
plicability of an equivalent current density in estimat¬ 
ing point-like and distributed current source configu¬ 
rations was investigated using simulated multichannel 
MCG data. 

In the literature, the calculation of the equivalent cur¬ 
rent density was first formulated in terms of the mini¬ 
mum norm estimate (MNE) [2], The problem of solv¬ 
ing the unknown source strengths leads to an ill-posed 
problem. The ill-posedness is a mathematical reflec¬ 
tion of physical phenomena that include the attenu¬ 
ation and smoothing effects of the volume conduc¬ 
tor, and the measurement noise. Rapid spatial vari¬ 
ations in the current distribution inside the heart are 
blurred and smoothed in the measured magnetic data. 
Therefore, a small amount of perturbation (noise, er¬ 
rors in the volume conductor model, discretization ef¬ 
fects etc.) in the MCG data tends to be magnified in 
the inverse solution. The calculation can be stabilized 
by using regularization. The regularization produces 
solutions which usually “fit” worse to the measured 
data than the original, unregularized solution but are 
more realistic and stable. 

2 Methods 

In the following, the calculation and the regulariza¬ 
tion methods used in this study to obtain current 


density images of the source current are presented. 
From the possible regularization schemes, zero- and 
second-order Tikhonov regularization were selected 
as means to obtain stable solutions. In addition, the 
benefits of using a pre-conditioning matrix in the cal¬ 
culation were evaluated. 

2.1 Lead fields 

The derivation of the equations is initiated by the con¬ 
cept of lead fields which can be thought as transfer 
functions between the current sources and the gener¬ 
ated magnetic signals. Therefore, the magnetic field 
values B t , i = !.•••. m, can be linked to the primary 
current distribution J p inside the heart via the lead 
fields L i 

Bi(vi) = J Lj(r ; ) ■ J p (r ; ) dr 1 , (1) 

where r , is the location of the 't h sensor. The inte¬ 
gration volume V' contains all current sources. Dis¬ 
cretization of Eq. 1 produces a matrix formulation 

B = LJ , (2) 

where vector B mx i contains the magnetic field 
values and matrix L mx 3 „ = (Li,---,L m ) is 
composed of the lead fields. Vector J. 3 .„ x i = 
(Ji iT . J] y . J\ z . J-ix■ ■ • •, J n z) T contains the unknown 
source strengths. The conventional minimum norm 
solution of an under-deterministic matrix equation, 
such as Eq. 2, can be expressed as 

J = Lt B = L T (LL T )- 1 B , (3) 

where L t is the pseudo-inverse of L. In practice, the 
unknown source strengths J cannot be directly solved 
from Eq. 3 because the lead fields of adjacent sensors 
are nearly linearly dependent on each other. There¬ 
fore, the inner product matrix {LL T ) has a low nu¬ 
merical rank, and the classical solution will be con¬ 
taminated by noise. The calculation can be stabilized 
by using regularization. 




2.2 Tikhonov regularization 

In T ik honov regularization, the form of the function 
to be minimized can be expressed as 

||LJ — B|| 2 q- A||ft!J|| 2 - Ymin , (4) 

where || • || is the Euclidean norm, R is the regular¬ 
ization operator, and A is a parameter which controls 
the weight given to the minimization of the side con¬ 
straint. Thereafter, the solution to the minimization 
problem in Eq. 4 can be written as 

J = L reg B = {L t L + \R t R) 'L t B , (5) 

or alternatively 

J = {R T R)- 1 L T (L(R T R)- l L T + AT) -1 B , (6) 

provided that the regularization operator R has a full 
rank 3 n. Both solutions involve the inversion of a 
3 n x 3 n matrix, therefore increasing the computa¬ 
tional demands with an increasing number of source 
points. 

2.2.1 Regularization operator 

The regularization operator R can be implemented in 
several ways. The use of the identity matrix I tends to 
minimize the norm of the solution. In that case, Eq. 6 
reduces to J = L 7 (/. L' f XI) 'B. which involves 
inversion of an m x m matrix only, and is therefore 
quickly computed. 

Discrete approximations of derivative operators (dif¬ 
ference, Laplacian) used as the regularization opera¬ 
tor tend to reduce either steepness or steep changes 
in the solution. The discrete surface Laplacian oper¬ 
ator A can be constructed, e.g., by using one of the 
approximations presented by Huiskamp [3]. 

2.2.2 Selecting a value for A 

In this study, the L-curve method was used to se¬ 
lect an optimal value for the regularization parameter 
A [4]. The name of the method relates to the shape 
of a continuous curve Co : {||LJ - B||, ||i2J||} re¬ 
lated to the minimization problem in Eq. 4. By us¬ 
ing the generalized singular value decomposition of 
the matrices L = UCX -1 and R = VSX -1 , the 
regularizing matrix L reg in Eq. 5 can be written as 
X[C T C + A S T S)~ 1 C T U T . In the decomposition, 
U and V are orthonormal, and X is non-singular. If 
L is assumed to have full rank m and R is assumed 


to have full rank 3n, the matrices C and S for the 
under-deterministic case in Eq. 2 are of the form 


C — ( O mx ( 3 „_ m ) C m ), 



where C m = diag(ai, • • ■, a m ) and S m = 
diag(/3i, ■ • •, p m ). The generalized singular values of 
matrices L and R are defined as m = ^. Therefore, 
the regularized solution J in Eq. 5 can be expressed 
as 

j = <7) 

where u, refers to the I th column of U and x, 3 „ rn+l 
to the 3 n—m+i th column of X. By using the expres¬ 
sion in Eq. 7, the squared norm of the residual can be 
written as 

||LJ-B|| 2 = f^-^( Ul -B)J , (8) 

and the squared norm of the side constraint is 

ll* J ll 2 = £[^K- B )] • ( 9 ) 

In the L-curve method, the goal is to find a value for 
A which balances the minimization between the resid¬ 
ual and the side constraint at the “comer” point of the 
curve. The comer is defined to have a maximal cur¬ 
vature. 

2.3 Weighted solutions 

In the calculation of weighted solutions, a pre¬ 
conditioning matrix was added into the original ma¬ 
trix formulation of Eq. 2: 

B = LDD 1 J . (10) 

If the product of matrices L and D is denoted as L, 
the solution of the under-determined Eq. 10 can be 
written as 

J = D[L r {LL r )- 1 B] . (11) 

In this study, the diagonal elements of the weighting 
matrix were selected to be the goodness of fit ( g ) val¬ 
ues for each source position. Despite the addition of 
the weighting matrix D, the matrix to be inverted in 
Eq. 11 is close to singular, and regularization needs to 
be used in obtaining a stable solution to Eq. 10. 
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Figure 1: Current density images of a shallow simu¬ 
lated current dipole, placed close to the anterior sur¬ 
face of the epicardium. The solutions were calculated 
with a) zero-order, b) second-order and c) weighted 
zero-order Tikhonov regularization. 


3 Results 

Simulated MCG data with a realistic amount of mea¬ 
surement noise were used to study different regu¬ 
larization operators. Point-like and distributed cur¬ 
rent sources were placed on the surface of or inside 
a realistic epicardium. A homogeneous BE torso 
model was used to model the volume conductor both 
in the forward and in the inverse calculations. The 
magnetic sensor configuration used in the simula¬ 
tions corresponded to the configuration of the 67- 
channel cardiomagnetometer at the BioMag Labora¬ 
tory at Helsinki University Central Hospital. 

3.1 Shallow and deep point sources 

In Fig. 1, current density imaging results of a shallow 
current dipole source are presented. The results were 
calculated with zero- and second-order T ik honov reg¬ 
ularization. The simulated dipole is situated 4 cm be¬ 
low the surface of the torso, approximately 10 cm un¬ 
der the plane of the sensors, close to the anterior sur¬ 
face of the right ventricle. The amplitude of the dipole 
moment was selected to be 10 pAm, and the direction 
was randomized. A realistic measurement noise with 
an RMS value of 80 fT was added to the simulated 
data. The resulting SNR was 15. From Fig. 1 it can 
be seen that all three calculation methods were able to 
localize a shallow dipolar source. The largest current 
density values were found close to the source, marked 
as a yellow arrow, on the anterior surface of the epi¬ 
cardium. 

In contrast to a shallow source, remarkable differ¬ 
ences were found for a deep current dipole source, 
situated 13 cm below the surface of the torso. This po¬ 
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Figure 2: Current density images of a deep simulated 
current dipole, placed close to the posterior surface 
of the epicardium. The solutions were calculated with 
a) zero-order, b) second-order and c) weighted zero- 
order Tikhonov regularization. 


sition was selected to represent a deep left ventricular 
source. The amplitude and the direction of the dipole 
moment were the same as for the shallow source. The 
current density images of the deep simulated source 
are presented in Fig. 2. It was found that only with 
second-order Tikhonov regularization, the source re¬ 
gion could be correctly localized. However, the re¬ 
gion with the largest current density amplitude was 
more wide than in the estimates in Fig. 1. The normal 
and the weighted zero-order regularization, in turn, 
resulted in false descriptions of the underlying source 
configuration. In these estimates, the largest current 
density values were found on the anterior surface of 
the epicardium. 

3.2 Distributed sources 

Distributed current sources were simulated by placing 
several current dipoles close to the anterior and to the 
posterior surface of the epicardium. The sources con¬ 
sisted of 16 current dipoles, each dipole having 1/16 
of the amplitude of the single dipole used in Figs. 1 
and 2. The dimension of the source area was approx¬ 
imately 4 cm. The current density images of the shal¬ 
low source region are presented in Fig. 3. The results 
closely resemble to the images in Fig. 1, however, the 
areas of large current density values are larger than for 
the single dipole source, corresponding well to the ex¬ 
tension of the source region. For a shallow distributed 
source, all the calculation methods were able to local¬ 
ize the current region correctly. 

The results for a deep distributed source are presented 
in Fig. 4. The phenomena found for a point source 
in Fig. 2 were reproduced in the current density im¬ 
ages of the deep distributed source. The normal zero- 
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Figure 3: Current density images of a shallow dis¬ 
tributed current source, consisting of 16 current 
dipoles spaced close to the anterior surface of the epi- 
cardium. The solutions were calculated with a) zero- 
order, b) second-order and c) weighted zero-order 
Tikhonov regularization. 
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Figure 4: Current density images of a deep dis¬ 
tributed current source, consisting of 16 current 
dipoles spaced close to the posterior surface of the 
epicardium. The solutions were calculated with a) 
zero-order, b) second-order and c) weighted zero- 
order Tikhonov regularization. 


order and weighted zero-order Tikhonov regulariza¬ 
tion resulted in false descriptions of the underlying 
source configuration where the largest current density 
values were found on the anterior surface of the epi¬ 
cardium. With second-order Tikhonov regularization, 
however, the deep source region could be correctly 
localized. In the second-order regularized image, the 
region with the largest current density amplitude cor¬ 
responded well to the true extension of the source re¬ 
gion. 

4 Discussion 

The results of this study showed that current den¬ 
sity imaging can localize both shallow and deep point 
sources but the regularization method applied in the 


calculation will have to be selected carefully. How¬ 
ever, the true power of the current density imag¬ 
ing technique is that it can also be used in localiz¬ 
ing several simultaneous point sources and distributed 
sources without the constraint of point-like source 
models, such as the current dipole. The results of 
the present study, as well as the previously obtained 
clinical results [5, 6, 7], encourage the use of current 
density imaging technique in cardiomagnetic source 
localization studies. 
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